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ABSTRACT 

A numerical  study  is  drne  comparing  three  algorithms  for  computing 
best  rational  approximations  using  the  uniform  (Chebyshev)  norm.  The 
experimental  results  and  theoretical  considerations  indicate  that  the 
Remes-dlfcor  algorithm  is  superior  as  a general-purpose  routine  to  both 
the  widely-used  Remes  algorithm  and  the  differential  correction  algorithm. 
The  three  algorithms  are  briefly  described  and  discussed,  and  the  experi- 
mental results  for  70  examples  are  presented  in  six  tables. 


1.  Introduction 

He  consider  the  problem  of  approximating  a given  real-valued  function 
f(x)  on  a finite  subset  T of  an  interval  [a,  b]  by  a rational  function 

R™<x>  • Pm(x)/Q  (x)  = [ I a,xJ]/[  ? b,xJ]  (1) 

n j«0  J j=0  J 

where  the  Integers  m and  n are  given,  (Pm,  Qn)  = 1,  Qn(x)  > 0 on  T,  and 

Q_(x)  Is  normalized  by  taking  max  lb. I = 1.  It  is  desired  to  choose 
n 0<j<n  J 

Rjj(x)  to  minimize  the  error  norm  ||f  - R™||  = max{|f(x)  - R™(x)|:  x eT}. 

We  note  that  if  f e C[a,  b]  and  the  best  approximation  to  f on  [a,  b]  Is 
nondegenerate  (i.e.  min(m  - aPm,  n - aQn ) » 0,  where  3P  * degree  of  P), 
then  the  best  approximation  on  T can  be  made  arbitrarily  close  everywhere 
on  [a,  b]  to  the  best  approximation  on  [a,  b]  by  taking  T sufficiently  dense 
in  [a,  b]  ([3]).  Thus  in  many  cases  the  continuous  problem  can  be  effectively 
treated  by  discretization  to  a finite  subset. 

Research  supported  in  part  by  National  Research  Council  of  Canada  under  grant 
A8061,  by  a University  of  Victoria  Faculty  Research  grant,  by  the  Air  Force 
Office  of  Scientific  Research,  Air  Force  Systems  Command,  USAF,. under  grant 
AF0SR-76-2878C,  and  by  the  National  Science  Foundation  under  grant  MCS-78-05847. 
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If  a computer  installation  wished  to  have  a single  library  routine 
for  calculating  best  uniform  rational  approximations,  then  according  to 
Lee  and  Roberts  [7]  they  should  select  either  the  Remes  or  differential 
correction  algorithm.  To  be  safe  they  should  select  the  differential 
correction  algorithm  because  of  its  guaranteed  convergence  properties; 
however,  the  costs  of  this  choice  are  a significantly  slower  routine  than 
the  Remes  algorithm  and  increased  storage  requirements.  The  main  purpose 
of  this  paper  is  to  provide  evidence  that  the  Remes-difcor  algorithm 
([4],  [8])  would  be  a better  choice  than  either  Remes  or  differential 
correction  as  a general  purpose  algorithm. 

The  three  algorithms  and  convergence  theorems  are  discussed  in 
Section  2 of  this  paper.  The  numerical  examples  are  given  in  Section  3, 
and  conclusions  are  drawn  in  Section  4. 


2.  The  Algorithms 

The  version  of  the  Remes  algorithm  considered  here  is  the  same  as 
considered  by  Lee  and  Roberts  ([7]).  This  algorithm  is  based  upon  alter- 
nation. At  the  completion  of  the  (k  - l)st  iteration  it  produces  a 
"reference  set"  Xk  * {xj[,  x^,  ....  xJj+n+2}  £ T-  In  the  kth  Iteration, 
the  algorithm  first  calculates  (If  possible)  the  best  approximation  Pk/Qk 
on  Xk  using  Newton's  method  to  solve  the  nonlinear  system  of  equations 
for  unknowns  aQ,  ...»  am,  b-j,  ....  bn  (the  coefficients  of  P and  Q)  and  x 
f(x*)  - P(xJ)/Q(x{)  * (-1 )1X,  1 « 1,  ....  m + n + 2,  with  bQ  = 1 (2) 

or,  equivalently, 

f(x})Q(x!{)  - P(x![)  = (-1  )ixQ(xJ),  1 -1,  .«,m  + n + 2,  with  bQ«l.  (3)  • 
If  Pk/Qk  is  not  the  best  approximation  on  T,  the  next  reference  set,  Xk+1, 
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Is  constructed  from  the  extreme  points  of  f - Pk/Qk.  This  is  done  by 
selecting  xj[+1  < ...  < xj^+2  each  in  T,  such  that  sgn(f (x^+1  J-P^Cx^1  )/Q|c(x|[+1 ) 

» -sgn(f(x£])  - Pk(x^])/Qk(xit]))»  1 = 1 + n + 1, 

.|f(xf+1)  - Pk(xJf+1)/Qk(x|f+1)|  > max|f(x{)  - pk(x- )/Qk(x- ) | , 1 = 1 

m + n + 2, and  for  some  n,  1 < n£m  + nt2,  |f(x[j+1)  - pk(xn+1 )/Qk(xn+1 )l 

* || f - Pk/Qk||  . This  is  normally  referred  to  as  a multiple  exchange. 

The  Initial  reference  set  is  chosen  to  be  the  points  of  T which  are  closest 
to  the  extreme  points  of  the  (m  + n + l)s^  Chebyshev  polynomial  translated 
to  [a,  b].  This  algorithm  can  fail  to  converge;  this  failure  can  be  caused 
by  the  desired  solution  having  a reference  set  of  less  than  m + n + 2 points, 
by  failure  to  be  able  to  solve  the  nonlinear  system  for  acceptable  Pk/Qk> 
or  by  making  an  improper  exchange  due  to  the  existence  of  poles  of  Pk/Qk 
off  of  the  reference  set  Xk>  It  has  been  shown  by  Burke  ([3]),  however, 
that  if  the  best  approximation  to  f on  [a,  b]  is  nondegenerate,  T is 
sufficiently  dense  in  [a,  b],  and  the  choice  of  initial  reference  set  is 
sufficiently  good,  then  the  algorithm  will  converge  to  the  best  approximation 
on  T. 

The  differential  correction  program  used  in  this  paper  is  an  improved 
version  of  a code  which  appeared  in  [6].  Recently,  another  version  of  [6] 
has  appeared  [2]  which  is  claimed  to  be  about  30%  faster  and  a little  more 
robust  than  the  code  In  [6].  In  order  to  describe  the  differential  correc- 
. tlon  algorithm,  let  us  assume  that  at  the  (k-l)st  step  an  approximation 
Pk/Qk  has  been  calculated  with  = ||f  - Pk/Qk||  . The  algorithm  then 
proceeds  to  calculate  the  next  approximation,  Pk+1/Qk+1,  by  using  linear 
programming  to  solve  the  problem 


L. 


I f (x)Q(x)  - P(x)|  - AkQ(x) 

minimize  max # . 

XeT  gku' 

.subject  to  |bj|  < 1,  j * 0,  1, ...»  n. 


The  Initial  approximation  Is  found  by  minimizing  max|fQ  - P|  with  bQ  = 1, 

xeT 

with  extra  constraints  added  to  ensure  that  Q > 0 on  T.  It  has  been  shown 
by  Barrodale,  Powell,  and  Roberts  ([1])  that  the  algorithm  has  guaranteed 
monotone  convergence;  that  is,  Ak  + inf  ||f  - R™||even  if  no  best  approximation 
on  T exists. 

The  Remes-difcor  algorithm  ([4],  [8])  is  a hybrid  of  the  two  algorithms 
described  above.  It  differs  from  the  Remes  algorithm  described  above  in 
two  crucial  respects.  First,  approximations  on  reference  sets  are  found 
using  the  differential  correction  algorithm  rather  than  by  solving  a nonlinear 
system  of  equations.  Thus,  an  approximation  with  a positive  denominator  on 
the  reference  set  is  guaranteed  even  if  the  system  has  no  solution.  Second, 

If  a g-pole  (that  is,  a point  where  the  denominator  is  very  small  in  absolute 
value  or  negative)  occurs  somewhere  in  T off  the  reference  set,  and  f - Pk/Qk 
changes  sign  m + n + 1 times  on  the  reference  set,  then  the  next  reference 
set  is  expanded  to  include  the  point  where  Qk  is  smallest.  Note  that  the 
flexibility  achieved  by  using  differential  correction  on  the  reference  set 
Is  essential  here,  since  the  new  reference  set  will  have  m + n + 3 points 
Instead  of  m + n + 2.  One  point  will  be  deleted  from  this  enlarged  reference 
set  after  Pk+]/Qk+]  is  computed.  The  Remes-difcor  algorithm  will  converge 
regardless  of  the  choice  of  initial  reference  set,  providing  that  g-pole 
free  best  approximations  exist  on  all  reference  sets  encountered. 

3.  Numerical  Examples 

The  three  algorithms  have  been  run  on  each  of  fourteen  functions  with 
given  sets  T (see  Table  I)  using  rational  functions  of  the  form  Pl^l  * 
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Pj/Q*.  P,  /Qj,  and  The  first  eleven  functions  are  the  functions  used 

by  Lee  and  Roberts  [7].  The  next  two  functions  were  designed  to  give 

113 

degenerate  approximations.  (Note  that  (0,  y,  U is  a set  of  alternating 

extreme  points  for  f12  (x)-l/(l  + x)  and  {0,  .1,  .2,  ...»  .9,  1}  is  a set  of 
alternating  extreme  points  for  f-jjU)  - 1 /[ (11/20)  - x]^.)  The  last  function 
provides  examples  where  best  approximations  do  not  exist. 

All  computations  were  done  with  the  University  of  Victoria  IBM  370/145 
computer  using  double  precision  arithmetic  (15  digits).  The  Remes  algorithm 
was  terminated  when  the  last  two  reference  sets  agreed.  The  differential 
correction  algorithm  was  terminated  when  (||f  - Pk/Q|J|  - ||f  - Pk+l^^k+1  H ^ 

|| f - Pfc/QiJI  was  <_  10  or  after  50  iterations.  The  Remes-difcor  algorithm 
was  terminated  when  there  were  no  g-poles  and  the  maximum  absolute  error 
on  the  reference  set  was  within  10’^  of  the  maximum  absolute  error  on 
the  reference  set,  or  after  20  exchanges. 

Table  II  contains  a sunmary  of  the  results.  In  Tables  Ilia,  b,  and  c 
we  give  the  CPU  times  (in  seconds),  and  point  out  any  examples  where  the 
computed  approximation  had  a pole  in  [a,  b].  In  Table  IV  we  give  the  error 
norms  computed  with  the  differential  correction  algorithm;  the  other 
algorithms  produced  the  same  error  norms  and  best  approximations  whenever 
they  produced  pole-free  approximations. 


4.  Conclusions 

From  Table  II  we  observe  that  Remes-difcor  was  nearly  as  robust  as 
differential  correction  and  much  more  robust  than  Remes  on  the  set  of 
examples  run.  The  one  failure  of  Remes-difcor  was  due  to  cycling  of  best 


approximations*  which  occurred  because  of  the  apparent  nonexistence  of  a 
best  approximation  on  some  of  the  reference  sets. 

Remes-difcor,  although  not  as  fast  as  Remes,  was  considerably  faster 
than  differential  correction.  Since  the  size  of  the  point  set  over  which 
approximations  are  actually  computed  at  each  stage  is  independent  of  the 
size  of  T for  Remes  and  Remes-difcor,  it  is  not  surprising  that  as  the  size 
of  T Increased  the  ratio  of  differential  correction  time  to  Remes-difcor 
time  Increased,  while  the  ratio  of  Remes-difcor  time  to  Remes  time  decreased. 

Differential  correction  does  have  the  advantage  of  flexibility;  since 
an  alternating  theory  is  not  required,  it  can  be  used  for  such  things  as 
simultaneous  approximation  (see  [5])  and  approximation  of  functions  of 
several  variables.  This  provides  ar.  additional  argument  for  Remes-difcor 
versus  Remes,  since  a Remes-difcor  program  contains  the  differential 
correction  subroutines,  which  can  be  called  directly  for  problems  on  which 
Remes-difcor  is  inapplicable  or  fails.  On  the  other  hand,  Remes-difcor 
can  handle  many  problems  for  which  differential  correction  cannot  be  used 
(because  of  storage  problems)  or  would  take  excessive  computer  time. 

In  summary,  Remes-difcor  combines  maximum  robustness  with  good  speed, 
and  appears  to  us  to  be  the  best  general  purpose  uniform  rational  approxi- 
mation algorithm  available  today.  A FORTRAN  listing  is  available  from  the 
first  author. 
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TABLE  I 

Data  Sets  Used  in  the  Numerical  Study 


Function  f(x) 

lVV 

Number  of  points 

1 

(Eoually  spaced) 

V ( 

X 

B 

[-1.11 

51 

sin(x) 

[-3,31 

21 

V 

/x 

10,11 

11 

1 

[0,0.5) 

V ' 

0 

x * 0.5 

21 

k 

-1 

(0.5,11 

V {« 

K 

0.5X+O.4 

[0,11 

(1,21 

51 

► 

f • < i 

• V .1 

[0,1] 

21 

D 

V 

e -e  +e 

(1,21 

V 

log(l+x) 

[0,11 

51 

ft  erf(x) 
o 

[0,2] 

21 

2 

-x 

a 

Hi 

•• 

[0,2] 

11 

£io:  r(x) 

[2,31 

51 

fll‘  ’ 

r(x) 

[2,3] 

101 

-24X/5+3/2 

[O.i] 

52x/15-17/30 

(4'2J 

£12!  1 

—92x/2 1+4 7/14 

<H> 

21 

26X/7-19/7 

k 

f 

(|.U 

(11/20- k/10)”2  + 0.1 (-l)k 

x - k/10 

£13S  1 

(k  - 0,1,. ..,10) 

11 

straight  line  joining 
adjacent  points 

x c [0,1],  x / k/10 

e 

10,1) 
x - 1 

21 

rl 


' 
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average  based  on  the  57  examples  in  which  all  three  algorithms  converged. 

average  based  on  the  23  examples  in  which  all  three  algorithms  converged 
and  T contained  51  or  more  points. 


Symbols  used  in  Tables  Ilia,  Illb,  IIIc 
F * algorithm  fails  to  converge. 

* * algorithm  converges  to  an  approximation  with  a pole  in  [a,  b], 
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TABLE  I 1. 1 a 
The  Remes  Algorithm 
CPU  Times  (in  seconds) 


P,/Q, 

VQ2 

P2/q2 

Pl/Q3 

V»2 

f1 

0.75 

0.72 

0.83 

0.82 

0.65 

'2 

0.24 

F 

0.43* 

0.50* 

0.30 

f3 

0.27 

0.17 

0.41 

0.39 

0.60 

f4 

0.20 

F 

0.56 

0.53 

0.92 

f5 

F 

0.98 

F 

1.21 

2.36 

U 

0.27 

0.40 

0.70 

0.58 

1.24 

f7 

0.47 

0.72 

0.59 

0.56 

1.04 

f8 

0.24 

0.24 

0.50 

0.46 

0.45 

U 

0.26 

0.20 

0.36 

0.39 

0.33 

f10 

0.73 

0.69 

0.56 

0.56 

0.76 

fll 

1.38 

1.36 

1.01 

1.56 

1.10 

f12 

0.58* 

F 

F 

0.60* 

F 

f13 

0.29 

0.38* 

1.16* 

F 

0.99* 

f14 

F 

F 

F 

F 

F 

TABLE 
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TABLE  I I I c 
The  Remes-difcor  Algorithm 
CPU  times  (in  seconds) 


V«i 

P0/Q2 

p2/q2 

V«3 

V«2 

*1 

0.86 

0.9/ 

1.82 

1.68 

2.43 

f2 

0.33 

0.10 

2.87 

6.42 

1.01 

f3 

0.72 

1.03 

1.93 

2.01 

4.51 

f4 

0.64 

0.18 

4.10 

2,53 

5.25 

f5 

. 1.06 

2.06 

4.66 

3.19 

13.31 

0.59 

0.81 

2.70 

3.95 

7.14 

f7 

0.46 

1.84 

1.13 

1.14 

3.39 

f8 

0.60 

0.87 

2.24 

2.10 

2.49 

f9 

0.79 

0.82 

2.33 

2.32 

2.95 

f10 

0.71 

0.90 

1.21 

3.30 

2.63 

f11 

0.78 

0.90 

1.30 

4.82 

2.26 

f12 

2.22 

1.21 

2.76 

2.97 

F 

f13 

0.68 

0.90* 

2.59* 

3.48* 

7.19* 

'l« 

1.08 

1.38 

4.01 

2.59 

4.95 

TABLE  IV 


Errors  of  Best  Approximation 


VO,  VQ2  V«2  VQ3  P4/<!2 


'l 

0. 20932 (-1) 

0. 34791 (-1) 

0.86644(-4) 

0.12392(-3) 

0.21 037 (-6 ) 

f2 

0.62542(0) 

0.99749(0) 

0.30608(0) 

0.30608(0) 

0.66482 (-2) 

f3 

0.36243H) 

0.18078(0) 

0.77019(-3) 

0.37281 (-2) 

0.10202(-4) 

f4 

0.81818(0) 

0.10000(1) 

0.26923(0) 

0.26923(0) 

0.70465(-l ) 

fS 

0.5891 6 (-1 ) 

0.22594(0) 

0.54260(-l ) 

0.48581 (-1) 

0.1 871 7 (-1 ) 

f6 

0.30872(0) 

0.20697(0) 

0.86504(-l ) 

0.95354(-l) 

0- 3091 9 (-1  ) 

f7 

0.85978(-3) 

0.92869(-l) 

0.17028(-5) 

0.74224 (-5) 

0.58255(-8) 

f8 

0.44085(-l ) 

0.19844(0) 

0. 1 3754(— 2 ) 

0.92931 (-3) 

0.44515(-4) 

f9 

0.72164(-1) 

0.69042(-l ) 

0.25586 (-2 ) 

0.41421 (-2) 

0.38213(-4) 

f10 

0.64376(-2) 

0.64307 (-2) 

0.36395 (-4) 

0. 55096 (-4 ) 

0.17428(-6) 

fn 

0.64420 (-2) 

0.64351 (-2) 

0.36432 (-4) 

0.551 60(— 4) 

0.1 7660(-6 ) 

f12 

0.50000(0) 

0.50000(0) 

0.50000(0) 

0.50000(0) 

0.10186(0) 

f13 

0.19754(+3) 

0.10000(0) 

0.10000(0) 

0.10000(0) 

0.10000(0) 

f14 

0.28934(-7) 

0. 31 263(-7) 

0.41851 (-7) 

0.1 6880 (—7 ) 

0.84017(-7) 
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